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Abstract 

By means of standard and histogram Monte Carlo simulations, we investigate the 
critical and compensation behaviour of a ternary mixed spin alloy of the type 
ABpCi^p on a cubic lattice. We focus on the case with the parameters correspond- 
ing to the Prussian blue analog (Ni^ Mn^Lj,)i.5[Cr^^^(CN)g] • nH20 and confront our 
findings with those obtained by some approximative approaches and the experi- 
ments. 
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1 Introduction 



Prussian blue analogs, a class of molecular-based magnetic materials, are of 
high interest due to a relatively easy control of their magnetic properties dur- 
ing a synthesis process through the selection of desired spin sources. These 
compounds exhibit some remarkable properties (see, e.g., [HE] and references 
therein), such as a photoinduced magnetic pole inversion, a compensation or 
multicompensation behaviour, or an anomalous thermal dependence of the 
susceptibility as the compensation temperature approaches the critical one. 
In particular, an example of the compensation behaviour is the Prussian blue 
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analog of the type (Ni"Mn"_p)i.5 [Cr™(CN)6] ■ nHaO, which includes both 
ferromagnetic {JNiCr > 0) and antiferromagnetic {JunCr < 0) superexchange 
interactions through the cyanide bridging ligands. This so called ternary metal 
compound, besides a critical temperature, exhibits a compensation tempera- 
ture for some range of the mixing ratio p of the metal ions, and an unusual 
behaviour of the susceptibility as the two temperatures merge [2l[3] . Thus, the 
compensation behaviour results from a proper selection of the parameter p, 
and its actual value also quantitatively affects the magnitudes of the critical 
and compensation temperatures. We note that the compensation behaviour, 
with important technological applications in the field of magneto-optic record- 
ing, is possible due to different temperature dependencies of the sublattice 
magnetizations. 

Due to their structural complexity, the theoretical investigations of the Prus- 
sian blue analogs have been mostly limited to a mean-field theory (MFT) 
treatment [3HZ] or an effective- field theory (EFT) approach [319], but also 
Monte Carlo (MC) simulations p^HT2] and exact recursion relations on the 
Bethe lattice [T3] have recently been employed. However, in the previous MC 
studies only simplified models of the Prussian blue analogs, consisting of a 
fixed spatial configuration of the spins placed on two equivalent (1:1 stoi- 
chiometry) interpenetrating sublattices, were considered. Furthermore, either 
lattice geometry ( [10]), spin values ( [12]), or both ( [TT]) did not correspond 
to the compound of our current interest. Therefore, the main goal of the 
present paper is to focus on a spin system corresponding to the Prussian blue 
analog (Ni^ Mn^'_p)i,5[Cr^^^(CN)g] ■ ?7,H20, which can be modeled by a three- 
dimensional ABpCi-p ternary alloy on a cubic lattice with a 3 : 2 stoichiome- 
try, consisting of three different Ising spins 5*^ = |, = 1, and Sc = ^ and 
the superexchange interaction ratio R = \Jac\/Jab = iJuncA/ Jmcr ~ 0.45 
[2111] • We also take into account fluctuations resulting from different spatial 
distributions of magnetic and non-magnetic sites. Further, we study how mag- 
netic properties of the ternary metal compound are modified when the param- 
eter R and the stoichiometry of the system are changed. 



2 Model and Monte Carlo simulations 



The mixed ferro-ferrimagnetic ternary alloy model of the type ABpCi-p is 
described by the Hamiltonian 



H=-Y: ^^Sf[JAB^,Sf + JAcSfil - 0)], (1) 
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where = ±|, ±1 for A ions, 5*^^ = ±1,0 for B ions, and 5"^^ = ±|, ±|, 
±1 for C ions, is a random variable which takes the value of unity (zero) if 
the site i is occupied by the A ion (vacant), and C,j is a random variable which 
takes the value of unity (zero) if the site j is occupied by the B (C) ion. To be 
consistent with the structure of the above-mentioned Prussian blue analog we 
also assume that the couplings between nearest neighbours include both the 
ferromagnetic {Jab > 0) and the antiferromagnetic {Jac < 0) interactions. A 
simulated Lx LxL lattice consists of two interpenetrating face-centered cubic 
sublattices, each one comprising sites. The A ions, which are randomly 
distributed on one sublattice with the concentration pa = 2/3, are alternately 
connected with the B or C ions, which are randomly located on the other 
sublattice with the concentration p or 1 — p, respectively. Thus, each A ion is 
coupled with six nearest-neighbouring ions of the type X {X = B or C) , while 
the ions of the type X have on average only four nearest neighbours of the type 
A. This arrangement reflects the 3:2 stoichiometry of the system. It should be 
stressed that in the previous MC studies (TOHIIKTI] the real structure of the 
Prussian blue analogs was not taken into account. 

We consider linear lattice sizes L ranging from 12 up to 50 with the periodic 
boundary conditions imposed. Initial spin states are randomly assigned and 
the updating follows the Metropolis dynamics. The lattice structure and the 
short range of the interactions enable vectorization of the algorithm. Since 
the spins on one sublattice interact only with the spins on the other, each 
sublattice can be updated simultaneously. Thus one sweep through the entire 
lattice involves just two sublattice updating steps. For thermal averaging, we 
consider up to 1 x 10^ MC sweeps in the standard and up to 1 x 10^ MC sweeps 
in the histogram MC simulations, discarding the first 20% of the total num- 
ber for thermalization. Finally, configurational averaging is performed over Z 
(typically 10 or 20) spatial configurations. Then the errors of the calculated 
quantities are determined from the values obtained for those configurations. 
Namely, the error bars in all the figures presented below will represent twice 
of the standard deviations. 



We calculate the sublattice magnetization per site 



rriA = [{Ma)] 



Na 
i=l 



(2) 



TUB = [{Mb)] 



L3' 



Nb 



T.sf 



(3) 



mc = [{Mc)] 



Nc 



c 



(4) 
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where Na is the number of A ions {Na = PaL^/'2) on one sublattice and Nb, 
Nq are the numbers of B and C ions {Nb = pL^/2, Nq = (1 —p)L^/2) on the 
other sublattice, and (■ ■ ■) denotes thermal and [■ ■ ■] configurational averages. 
The total magnetization per site is defined as 



m = [(M)] 



-(Ma + Mb + Mc 



1 



{niA + mB + mc) 



(5) 



and the staggered magnetization per site is given by 



(6) 



The corresponding direct and staggered susceptibilities per site xm and xMs 
respectively, are given by 



Xm 



((M2) - (M)2) 



(7) 



and 



knT 



{{Ml) - (M.)2) 



(8) 



Finally, the specific heat per site c is defined by 
\{H^) - {Hf) 



c — 



L^kBT^ 



(9) 



In the histogram MC simulations, developed by Ferrenberg and Swendsen p31 
[T6] , we perform longer runs at a chosen temperature and calculate the energy 
histogram P{H), the order parameters histograms P{0) {O = M or Ms), 
as well as the physical quantities ([2]) — ([9]). Using data from the histograms 
one can calculate physical quantities at neighbouring temperatures, and thus 
determine the values of extrema of various quantities and their locations with 
high precision for each lattice size. In such a way we can obtain quality data 
for finite-size scaling (FSS) analysis. The extrema of the direct and staggered 
susceptibilities at a second-order transition are known to scale with the lattice 
size as 

XO,maAL) OC L^o/uo ^ (10) 



where uq, 7o represent the correlation length and susceptibility critical expo- 
nents, respectively, pertaining to the order parameter O (O = M or Mg). 
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3 Results and discussion 

3.1 Critical temperature 

The critical temperature of the current model can be considerably affected 
by the values of the exchange ratio R and the concentration p. Phase dia- 
grams in a wide range of the (p — R) parameter space are estimated based 
on monitoring of the positions of the (critical) maxima in different thermo- 
dynamic quantities, such as xm, XMsi and c. Due to the extensive coverage it 
would be computationally too demanding to perform the FSS analysis in each 
considered point of the {p — R) plane to determine the critical temperatures 
with high precision. Therefore, we estimate the values of from the peaks 
positions of xm^ fo^' a fixed value of L, i.e., we determine the pseudocritical 
temperature Tc{L) . After checking in a few selected points we chose the value 
of L = 20 above which the values of Tc{L) do not tend to change significantly. 
In Fig. [H we show dependence of such estimated critical temperature on R 
for some selected values of p. 

[Fig. 1 about here.] 

The limiting cases of p = and p = 1 correspond to the binary mixed systems 
AC with spins | and |, and AB with spins | and 1, respectively. While, the 
latter binary system is not affected by the parameter R and gives a fiat line, 
the former one gives the transition line with the steepest gradient. Further, 
from the figure it can be seen that all the curves appear to intersect in a 
single point corresponding to the critical value Re- By closer inspection of 
the intersection area for various lattice sizes we locate Rc ~ 0.47, a value 
surprisingly close to that obtained by MFT, Rf^^ = (^)^/^ ^ 0.478 [6]. 
This finding comes surprising particularly in the light of the fact that the 
jlMFT yq\^q does not depend on the coordination numbers (numbers of nearest 
neighbours) in the respective sublattices, while the current approach accounts 
for the system's 3:2 stoichiometry. This prompted us to check the value of 
Rc in our system with different stoichiometries. In the considered cases of 
1 : 1 {pa = 1), 3 : 2 {pA_ = 2/3), and 3 : 1 {pa = 1/3) stoichiometries 
we find that the critical temperatures decrease with decreasing pa (hence, 
decreasing overall coordination number), however, the respective values of 
Rc are indistinguishable within the current precision. Moreover, they are very 
close to the value obtained by the MC simulations on a two-dimensional system 
[TU] . These results indicate very little or no sensitivity of the value of Rc to 
the lattice stoichiometry, supporting the MFT results. However, the value 
of the corresponding critical temperature obtained by the MC approach is 
considerably lower than that obtained by MFT. For example, in the system 
with the 3 : 2 stoichiometry the value of the critical temperature from the MC 
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calculation is ksT^^/JAB ~ 3.57, while MFT gives kBT^^^/JAB ~ 4.47. 

The little sensitivity of the value of Rc to the lattice stoichiometry can be 
understood considering the processes that the respective parameters control. 
The parameter R = \Jac\/ Jab controls the relative strength of the exchange 
interactions between A ions and their neighbours on C and B sublattices. At 
Rc, a kind of a balance sets in at which a mutual exchange of B and C ions 
does not affect the critical temperature of the system. On the other hand, 
the lattice stoichiometry (controlled by the parameter p^) relates solely to the 
concentration of A ions. Therefore, while the stoichiometry affects the overall 
coordination number and hence the magnitude of the critical temperature, it 
is not likely to affect the relative strength of the exchange interactions between 
the AB and AC sublattice pairs, i.e. the value of Rc- 

As we can observe from Fig. [H the transition temperature for the ternary 
systems with R < Rc is always lower than that of the binary mixture AB and 
increases with increasing concentration p, while the opposite trend is observed 
for R > Rc- This behaviour of Tc is evident from Fig. [2l 

[Fig. 2 about here.] 

in which the transition temperature is shown as a function of a varying concen- 
tration. The transition line with R = 0.45 corresponds to the real compound 
(NijMn"_p)i.5 [Cr"^(CN)6]-nH20 0] and, as the Fig. shows, due to its proxim- 
ity to the value of Rc, the critical temperature is little affected by the change 
of the mixing ratio p. The transition curves merge at p = 1, for which the 
ternary system becomes the binary mixture AB, and therefore, R becomes a 
redundant variable. 

An interesting situation arises at the low-concentration side for the vanishing 
value of R. Namely, the transition temperature decreases with the decrease 
of R but as long as R remains finite the phase transition occurs in the whole 
concentration range. The situation for R = corresponds to configurations 
involving "loose" ions of the types C and A with no mutual exchange interac- 
tions and their concentration increases with decreasing value of p. Therefore, 
even though we do not deal with a problem of a dilution with non-magnetic 
impurities, the system involves a portion of spins that in the absence of the 
exchange interactions behave like non-magnetic atoms, since they are inactive 
in producing magnetic ordering. For this reason we can expect some "percola- 
tion threshold" Pc below which no phase transition can occur. Unfortunately, 
within the current approach it is computationally very demanding to locate 
the value of pc with a higher precision. The main reasons are extremely long 
thermalization times near pc, necessity to consider large lattice sizes, and to 
average over a large number of configurations. For instance, in our tests on 
individual configurations for p = 0.2, the equilibrium was not reached even 
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after ^ 10^ MC sweeps. Nevertheless, from the observation of some physical 
quantities, as well as the finite-size scaling, it appears that the long-range 
ordering survives down to at least p = 0.3. In Fig. Owe show temperature de- 
pendence of the specific heat for p = 0.5, 0.3 and 0.1, for the maximum linear 
size of L = 32, in order to suppress as much as possible the finite-size effects. 
The specific heat was chosen due to the fact that it exhibits much smaller 
fluctuations, and hence "smoother" behaviour, than the susceptibility. While 
the curves for p = 0.5 and 0.3 in Fig. |3] display well-defined sharp maxima, 
typical for a second-order phase transition, only a broad maximum is observed 
for p = 0.1. 

[Fig. 3 about here.] 
[Fig. 4 about here.] 

To see qualitative difference between the anomalies in the physical quantities 
for p = 0.3 and 0.1, we additionally performed FSS analysis of the direct 
susceptibility xm (for R = 0, there is no reason to consider XmJ using the 
histogram MC simulation. As depicted in Fig.HJ for p = 0.3 the standard Ising 
value of the critical exponents ratio 7/z/ 1.97 [17] was recovered with a fair 
precision, while for p = 0.1 even a negative value was obtained. We conjecture 
that the anomalies in the response functions for p = 0.1 originate from the 
presence of finite clusters of ordered spins rather than the long-range ordering, 
and therefore, we roughly estimate the value of Pc to fall into the region of 

0. 1 < Pc < 0.3. On the other hand, MFT is not able to deal correctly with the 
percolation problem and predicts Pc = [6]. 

3.2 Compensation temperature 

[Fig. 5 about here.] 
[Fig. 6 about here.] 

The compensation temperature is the temperature at which total magne- 
tization of the system m vanishes below its critical temperature and it is 
determined by the condition m{Tk) = 0. In this subsection we do not inves- 
tigate in detail the compensation behaviour of the current model for various 
values of the parameters R and p but rather focus on the situation correspond- 
ing to the real Prussian blue analog (Ni^ Mn^^_p)i.5 [Cr"^(CN)6] ■ nB.20 case, 

1. e. the case of R = 0.45, for which we can confront our results with some 
experimental findings [3lll]. For this value of the exchange ratio we found a 
compensation behaviour for the concentration values p ranging approximately 
between 0.305 and 0.425, which roughly coincides with the results predicted 
by MFT [12] and EFT [H]. Some examples of the variation of the total mag- 
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netization with temperature in this concentration region are shown in Fig. |5l 
The points at which the magnetization curves change sing below the respective 
critical temperatures represent the compensation temperatures. The variation 
of the critical (filled circles) and compensation (empty circles) temperatures 
with the concentration p is presented in Fig. [61 These theoretical results are 
in line with the experimental magnetization measurements jl], in which the 
compensation behaviour was observed for the concentrations 0.38 < p < 0.42. 
The lower bound of p was not explored in this experimental study, but the 
susceptibility measurements in [3] suggest that the concentration at which the 
compensation line joins the critical one is close to p = 0.30 i.e., close to the 
theoretical predictions. 

[Fig. 7 about here.] 
[Fig. 8 about here.] 

The latter measurements revealed an unusual behaviour of the susceptibility 
at a critical temperature in the vicinity of the compensation point at p = 
0.30. Namely, for p = 0.30 the paramagnetic susceptibility did not diverge 
near T^. Motivated by this finding, we simulated the susceptibility for the 
same concentration and compared its behaviour with that further away from 
p = 0.30. The results are shown in Fig. [71 While the curves for p = 0.1 
and 0.5 show prominent peaks at the transition (proper divergence cannot 
be observed due to finite-size effects), the one for p = 0.3 is almost flat. 
We note that a signiflcant contribution to the relatively large error bars at 
the susceptibility peaks for p = 0.1 and 0.5, besides the fluctuations of the 
peaks' heights, come also from the fluctuations of the peaks' positions obtained 
from different realizations. On the other hand, the principal source of the 
errors away from the peaks, where normally relatively small fluctuations are 
expected, is decreasing reliability of the reweighted values from the HMC 
method as we move away from the simulation temperature (usually chosen 
close to the peak's location [151 [16]). Nevertheless, the small error bars for 
p = 0.3 in the vicinity of the critical temperature indicate absence of the 
characteristic peak in all realizations. 

This phenomenon can be explained by observing temperature dependences of 
the magnetization curves in this region of the concentration p. As shown in 
Fig. [5l for p > 0.3 the compensation effect occurs due to different tempera- 
ture dependences of the sublattice magnetizations. Namely, as the temperature 
approaches the critical value from below, the positive sublattice magneti- 
zations {rriA, ^b) prevail over the negative one (mc), resulting in the flnite 
(positive) value of the total magnetization. However, as p approaches the value 
of 0.3 (from either side), the positive and negative sublattice magnetizations 
tend to cancel out just below Tc (see the case oi p = 0.3 in Fig. [5]). This 
situation is characteristic for an antiferromagnetic system which exhibits zero 
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total magnetization and non- divergent behaviour of the direct susceptibihty 
at critical (or Neel) temperature. 

In order to explore this behaviour in more detail, we focused on the behaviour 
of the susceptibility critical exponents in the area where the critical and com- 
pensation temperature lines merge. The exponents are calculated from the 
finite-size scaling given by equation ( ITOj) . using the histogram MC simula- 
tions. Fig. [8] demonstrates that while the ratio of the exponents ^m^I'^Msi 
corresponding to the staggered susceptibihty, retains the values close to the 
standard Ising value 7/z/ ^ 1.97 [T7] in the whole concentration range, the 
ratio ^mI'^m-, corresponding to the direct susceptibility, clearly deviates from 
this value in a relatively broad range of p. Namely, it gradually decreases as 
the value of p = 0.3 is approached from either side and it practically van- 
ishes at p ~ 0.3. Thus, this MC result clearly confirms the experimental [3] 
and MFT [TS] predictions that the direct susceptibility of the Prussian blue 
analog (Ni^ Mn5^'Lp)i.5[Cr^"(CN)g] ■ nH20 for a special value of for which 
the compensation line terminates at the phase transition boundary, does not 
diverge at the critical temperature. Furthermore, it suggests that the suscep- 
tibility critical exponents of the model with a compensation behaviour might 
be affected by the proximity of the compensation point. This will be subject 
to further investigation. 



4 Conclusions 

We employed standard and histogram Monte Carlo simulations to study the 
critical and compensation properties of a ternary ferro-ferrimagnetic spin alloy 
of the type ABpCx^^ on a cubic lattice, focusing on the case with the lattice 
stoichiometry and the parameters corresponding to the Prussian blue analog 
(NiJ,^Mn^/_p)i.5[Cr"^(CN)g] ■ nHaO. We confirmed the existence of the critical 
superexchange interaction ratio Rc for which the critical temperature does not 
depend on the mixing ratio of the ions on B and C sublattices, which had been 
predicted within the MFA [6J as well as MC simulations on a square lattice 
with equivalent sublattices (1 : 1 stoichiometry) [ID]. More importantly, we 
found that the value of i?c ~ 0.47 is virtually insensitive to the change in 
the stoichiometry and, therefore the lattice coordination number, and is close 
to that predicted by MFT on the same lattice |E] and that evaluated by MC 
simulations on the square lattice [10]. In contrast to these results, a recent 
MC study on the present spin system with equivalent sublattices gave quite 
different value of Rc = 0.513 [IS]. However, we believe that this result is 
erroneous. Our suspicion is based not only on the deviating value of Rc, but 
also an unusually high value of the critical temperature for Rc or equivalently 
for p = 1. The latter corresponds to the binary mixture of the spins S'a = | 
and S'b = 1, and its critical temperature had already been established within 
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MFT [20] and EFT with correlations [21] to be respectively roughly equal 
and lower than that obtained by the MC simulation in [19]. However, it is 
well known that effective field theories generally overestimate the value of the 
critical temperature in comparison with MC simulations. 

Finally, we showed that for R = the long-range order does not survive down 
to p = 0, as predicted by MFT, and we roughly estimated the bounds for the 
percolation threshold pc- Focusing on the case of i? = 0.45 and 3 : 2 stoichiom- 
etry, corresponding to the real compound (Ni^ Mn^'Lp)i.5[Cr™(CN)g] ■ nH20, 
we reproduced the experimentally observed phenomenon of a non-divergent di- 
rect susceptibility at the critical temperature in the vicinity of a compensation 
point. This peculiar behaviour of the direct susceptibility was demonstrated 
by evaluation of the critical exponents in the finite-size scaling analysis. 

Based on the above considerations, it can be concluded that the current Monte 
Carlo study brought both quantitative and qualitative improvement over the 
results previously obtained within MFT. Nevertheless, despite its simplicity, 
MFT in many respects did a good job of predicting the critical and com- 
pensation behaviour of the ternary alloy considered in the study. This finding 
corroborates that resulting from the experimental study of the magnetic prop- 
erties of the Prussian blue analogs [2]. 
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Fig. 3. Temperature dependences of the specific heat for i? = 0, L = 32, and different 
values of p. 
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Fig. 5. Temperature dependences of the total magnetization m ioT R = 0.45 and 
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Fig. 6. Critical (•) and compensation (o) temperatures as functions of p for R = 0.45. 
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Fig. 8. Ratios of the critical exponents ^m^/^Ms (bullets) and ^mI^m (diamonds), 
obtained from finite-size scaling analyses of the staggered and direct susceptibilities, 
respectively, as functions of the mixing ratio p for i? = 0.45. 
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